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Abstract 

(N 

We show that the sensor self-localization problem can be cast as a static parameter estimation 
problem for Hidden Markov Models and we implement fully decentralized versions of the Recursive 
Maximum Likelihood and on-line Expectation-Maximization algorithms to localize the sensor network 
^\ , simultaneously with target tracking. For linear Gaussian models, our algorithms can be implemented 

exactly using a distributed version of the Kalman filter and a novel message passing algorithm. The 
latter allows each node to compute the local derivatives of the likelihood or the sufficient statistics 
• needed for Expectation-Maximization. In the non-linear case, a solution based on local linearization in 

the spirit of the Extended Kalman Filter is proposed. In numerical examples we demonstrate that the 
^\ • developed algorithms are able to learn the localization parameters. 

Collaborative tracking, sensor localization, target tracking, maximum likelihood, sensor networks 

1 Introduction 

This paper is concerned with sensor networks that are deployed to perform target tracking. A network 
■ is comprised of synchronous sensor-trackers where each node in the network has the processing ability to 
CN ! perform the computations needed for target tracking. A moving target will be simultaneously observed 
by more than one sensor. If the target is within the field-of-view of a sensor, then that sensor will collect 
measurements of the target. Traditionally in tracking a centralized architecture is used whereby all the 
. sensors transmit their measurements to a central fusion node, which then combines them and computes the 
t— ( | estimate of the target's trajectory. However, here we are interested in performing collaborative tracking, but 
without the need for a central fusion node. Loosely speaking, we are interested in developing distributed 
tracking algorithms for networks whose nodes collaborate by exchanging appropriate messages between 
neighboring nodes to achieve the same effect as they would by communicating with a central fusion node. 

A necessary condition for distributed collaborative tracking is that each node is able to accurately 
determine the position of its neighboring nodes in its local frame of reference. (More details in Section [2j) 
This is essentially an instance of the self -localization problem. In this work we solve the self-localization 
problem in an on-line manner. By on-line we mean that self-localization is performed on-the-fly as the 
nodes collect measurements of the moving target. In addition, given the absence of a central fusion node 
collaborative tracking and self-localization have to be performed in a fully decentralized manner, which 
makes necessary the use of message passing between neighboring nodes. 

There is a sizable literature on the self-localization problem. The topic has been independently pursued 
by researchers working in different application areas, most notably wireless communications (TJ [21 [3l EJ E] • 
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Although all these works tend to be targeted for the application at hand and differ in implementation 
specifics, they may however be broadly summarized into two categories. Firstly, there are works that rely 
on direct measurements of distances between neighboring nodes [2j |3j HI E] . The latter is usually estimated 
from the Received Signal Strength (RSS) when each node is equipped with a wireless transceiver. Given 
such measurements, it is then possible to solve for the geometry of the sensor network but with ambiguities 
in translation and rotation of the entire network remaining. These ambiguities can be removed if the 
absolute position of certain nodes, referred to as anchor nodes, are known. Another approach to self- 
localization utilizes beacon nodes which have either been manually placed at precise locations, or their 
locations are known using a Global Positioning System (GPS). The un-localized nodes will use the signal 
broadcast by these beacon nodes to self- localize [TJ [6j [7j [8]. We emphasize that in the aforementioned 
papers self- localization is performed off-line. The exception is [Sj, where they authors use Maximum 
Likelihood (ML) and Sequential Monte Carlo (SMC) in a centralized manner. 

In this paper we aim to solve the localization problem without the need of a GPS or direct measurements 
of the distance between neighboring nodes. The method we propose is significantly different. Initially, the 
nodes do not know the relative locations of other nodes, so they can only behave as independent trackers. 
As the tracking task is performed on objects that traverse the field of view of the sensors, information is 
shared between nodes in a way that allows them to self-localize. Even though the target's true trajectory 
is not known to the sensors, localization can be achieved in this manner because the same target is being 
simultaneously measured by the sensors. This simple fact, which with the exception of [9J [TUJ, [TT] seems 
to have been overlooked in the localization literature, is the basis of our solution^. However, our work 
differs from [10] in the application studied as well as the inference scheme. Both [U [10] formulate the 
localization as a Bayesian inference problem and approximate the posterior distributions of interest with 
Gaussians. |10| uses a moment matching method and appears to be centralized in nature. The method in 
[9] uses instead linearization, is distributed and on-line, but its implementation relies on communication 
via a junction tree (see |13| for details) and requires an anchor node as pointed out in [141 Section 6.2.3]. 
In this paper we formulate the sensor localization problem as a static parameter estimation problem for 
Hidden Markov Models (HMMs) |15[ PTE] and we estimate these static parameters using a ML approach, 
which has not been previously developed for the self-localization problem. We implement fully decentralized 
versions of the two most common on-line ML inference techniques, namely Recursive Maximum Likelihood 
(RML) [HllTUQjj] and on-line Expectation-Maximization (EM) [2QI EIJ [22]. A clear advantage of this 
approach compared to previous alternatives is that it makes an on-line implementation feasible. Finally, 

is based on the principle shared by our approach and [H [10] . In the authors exploit the correlation 
of the measurements made by the various sensors of a hidden spatial process to perform self-localization. 
However for reasons concerned with the applications being addressed, which is not distributed target 
tracking, their method is not on-line and is centralized in nature. 

In the signal processing literature for sensor networks one may find various related problems. In [23j 
a distributed EM algorithm was developed to estimate the parameters of a Gaussian mixture used to 
model the measurements of a sensor network deployed for environmental monitoring (see |24| for an on- 
line version.) In |25| a similar problem is treated using a distributed gradient method. We emphasize 
that in each of these papers the measurements correspond to a static source instead of a dynamically 
evolving target. In addition, a related problem is that of sensor registration, which aims to compensate for 
systematic biases in the sensors and has been studied by the target tracking community |26U27| . However, 
the algorithms devised in [261 [27] are centralized. Yet another related problem is the problem of average 
consensus |28| . The value of a global static parameter is measured at each node via a linear Gaussian 
observation model and the aim is to obtain a maximum likelihood estimate in a distributed fashion. Note 
that all the aforementioned papers, except [9] and [10] , do not deal with a distributed localization and 
tracking task. 

The structure of the paper is as follows. We begin with the specification of the statistical model for the 
A short preliminary version of the this work was published in the conference proceedings |12| . 
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localization and tracking problem in Section [2j In Section |3] we show how message passing may be utilized 
to perform distributed filtering. In Section U we derive the distributed RML and on-line EM algorithms. 
Section [5] presents several numerical examples on small and medium sized networks. In Sections [6] we 
provide a discussion and a few concluding remarks. The Appendix contains more detailed derivations of 
the distributed versions of RML and EM. 

2 Problem Formulation 

We consider the sensor network (V, £) where V denotes the set of nodes of the network and £ is the set 
of edges (or communication links between nodes.) We will assume that the sensor network is connected, 
i.e. for any pair of nodes i,j 6 V there is at least one path from i to j. Nodes i,j G V are adjacent 
or neighbors provided the edge (i,j) £ £ exists. Also, we will assume that if € £, then € £ 

as well. This implies is that communication between nodes is bidirectional. The nodes observe the same 
physical target at discrete time intervals n £ N. We will assume that all sensor-trackers are synchronized 
with a common clock and that the edges joining the different nodes in the network correspond to reliable 
communication links. These links define a neighborhood structure for each node and we will also assume 
that each sensor can only communicate with its neighboring nodes. 

The hidden state, as is standard in target tracking, is defined to comprise of the position and velocity 
of the target, X r n = LY^(l), X^(2), X^(3), A^(4)] T , where X^(l) and X^(3) is the target's x and y position 
while A^(2) and A^(4) is the velocity in the x and y direction. Subscript n denotes time while superscript 
r denotes the coordinate system w.r.t. which these quantities are defined. For generality we assume that 
each node maintains a local coordinate system (or frame of reference) and regards itself as the origin (or 
center of) its coordinate system. 

As a specific example, consider the following linear Gaussian model: 

X r n = A n X r n _ 1 + b r n + V n , n>l, (1) 

where V n is zero mean Gaussian additive noise with variance Q n and b r n are deterministic inputs. The 
measurement made by node r is also defined relative to the local coordinate system at node r. For a 
linear Gaussian observation model the measurement is generated as follows: 

y; = « + <f n + c, n>l, (2) 

where is zero mean Gaussian additive noise with variance R r n and d r n is deterministic. Note that the 
time varying observation model {(C r n , d r n , i?,^)} n >i is different for each node. A time-varying state and 
observation model is retained for an Extended Kalman Filter (EKF) implementation in the non-linear 
setting to be defined below. It is in this setting that the need for sequences {6^} n >i an d {d r n } n >i arises. 
Also, the dimension of the observation vector Y£ need not be the same for different nodes since each node 
may be equipped with a different sensor type. For example, node r may obtain measurements of the 
target's position while node v measures bearing. Alternatively, the state-space model in {J])-® can be 
expressed in the form of a Hidden Markov Model (HMM): 

KIK-l =<-l~/n(.K-l), (3) 

y;|x; = <~ 5 ;(.K), (4) 

where f n denotes the transition density of the target and g r n the density of the likelihood of the observations 
at each node r. 

Figure [T] (a) illustrates a three node setting where a target is being jointly observed and tracked by 
three sensors. (Only the position of the target is shown.) At node 1, X\ is defined relative to the local 
coordinate system of node 1 which regards itself as the origin. Similarly for nodes 2 and 3. We define 
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(a) Three node joint tracking example (b) Joint tracking error vs number of nodes 



Figure 1: Left: a three node network tracking a target traversing its field of view. The trajectory of the 
target is shown with the solid line. Each node regards itself as the center of its local coordinate system. 
At time n a measurement is registered by all three nodes. The ellipses show the support of the observation 
densities for the three nodes, i.e. the support of <?*(Y^|.) is defined as all x\ such that g^Y^x^) > ; 
similarly for the rest. The filtering update step at node 1 will clearly benefit from the observations made by 
nodes 2 and 3. The localization parameters 6l' 2 , 0*' 3 are the coordinates of node 1 in the local coordinate 
systems of node 2 and 3 respectively. While was defined to be the state of the target, which includes 
its velocity, for this illustration only, is to be understood as the position of the target at time n w.r.t. 
the coordinate system of node r. Right: Average absolute tracking error is plotted against the number 
of nodes to illustrate the benefit of collaborative tracking. The results are obtained using a centralized 
implementation with 50 independent runs, 10 4 time steps for a chain sensor network of different length 
and A n = B n = Q n = & n = D\ = R n = 1, b\ = < = 0. 

Q l * to be the position of node i in the local coordinate system of node j. This means that the vector X % n 
relates to the local coordinate system of node j as follows (see Figure [T]) : 

xi = xi + ei'i. 

The localization parameters are static as the nodes are not mobile. We note the following 

obvious but important relationship: if nodes i and j are connected through intermediate nodes j\ , ji , ■ ■ ■ , jm 
then 

0*J = 0jii + + 0jw'a + . . . + ei m - ujrn + 9i m ' j . (5) 

This relationship is exploited to derive the distributed filtering and localization algorithms in the next 
section. We define so that the dimensions are the same as the target state vector. When the state 
vector is comprised of the position and velocity of the target, only the first and third components of Ol' 1 
are relevant while the other two are redundant and set to 6>l J (2) = and 6^(4) = 0. Let 

0* = ^ = 0, (6) 

where 9l' z for all i G V is defined to be the zero vector. 



4 



Let Y n denote all the measurements received by the network at time n, i.e. Y n = {Yn}v£V- 

We 

also denote the sequence (Y%, ...,Y n ) by Y\- n . In the collaborative or joint filtering problem, each node r 
computes the local filtering density: 

PeS X n\ Y l:n) « P0„ Q^K)pS, «|*l:n-l), (7) 

where p r g^{x r n \Yi- n _i) is the predicted density and is related to the filtering density of the previous time 
through the following prediction step: 

P r eS X n\ Y l.n-l) = j /n«|<-l)P^«_l|>l : n-l)rf<-l- (8) 

The likelihood term is 

pl(YnK) = U9n( Y nK + 0l' V ), (9) 

veV 

where the superscript on the densities indicate the coordinate system they are defined w.r.t. (and the 
node the density belongs to) while the subscript makes explicit the dependence on the localization param- 
eters. Let also Mn|n~l and /i^ denote the predicted and filtered mean of the densities pg^ (a;^|Yi :n _i) and 
Pg (x r n \Yi :n ) respectively, where the dependence on 9* is suppressed in the notation. The prediction step in 
([8]) can be implemented locally at each node without exchange of information, but the update step in (|7|) 
incorporates all the measurements of the network. Figure Q] (a) shows the support of the three observation 
densities as ellipses where the support of g\iY^\-) is defined to be all x 1 such that g\iY^\ ) > 0; similarly 
for the rest. The filtering update step at node 1 can only include the observations made by nodes 2 

12 13 

and 3 provided the localization parameters 6* and 9*' are known locally to node 1, since the likelihood 
Pg (Y n \Xn) defined in ([9]) is 

Sn( Y n \ x n)9n( Y n \ x n @* )9n(Yn \ x n @* )■ 

The term joint filtering is used since each sensor benefits from the observation made by all the other 
sensors. An illustration of the benefit w.r.t. the tracking error is in Figure Q] (b). We will show in Section 
[3] that it is possible to implement joint filtering in a truly distributed manner, i.e. each node executes 
a message passing algorithm (with communication limited only to neighboring nodes) that is scalable 
with the size of the network. However joint filtering hinges on knowledge of the localization parameters 
0* which are unknown a priori. In Section H] we will propose distributed estimation algorithms to learn 
the localization parameters, which refine the parameter estimates as new data arrive. These proposed 
algorithms in this context are to the best of our knowledge novel. 

2.1 Non-linear Model 

Most tracking problems of practical interest are essentially non-linear non-Gaussian filtering problems. 
SMC methods, also known as Particle Filters, provide very good approximations to the filtering densities 
|29| . While it is possible to develop SMC methods for the problem presented here, the resulting algorithms 
require significantly higher computational cost. We refer the interested reader to |141 Chapter 9] for more 
details. In the interest of execution speed and simplicity, we employ the linearization procedure of the 
Extended Kalman filter (EKF) when dealing with a non-linear system. Specifically, let the distributed 
tracking system be given by the following model: 

K = MK-i) + v n , (io) 
K = 1MK) + wZ, (ii) 

where <f> n : M 4 — > M 4 and i\) r n : M 4 — > M. dy are smooth continuous functions. At time n, each node will 
linearize its state and observation model about the filtered and predicted mean respectively. Specifically, 



5 



Algorithm 1 Generic message passing at time n 



1: begin 

2: At k = 1, compute: 



<i = (14) 
= /-X'. (15) 



3: for k = 2, K compute: 



'■.7 

n,k 



F l + 



m 



p. I 

n,k—V 



~h3 



m 



n,k 



p£ne(i)\{j} 

pene(i)\{j} 



n,fc— 1' 



(16) 
(17) 



4: endfor 
5: end 



a given node r will implement: 

K = Mv r n-i) + vMv r n -i)(K-i - + K, (12) 

^ = C^ln-l) + VP n <ti\ n -l)(K - /£|„-l) + Wl (13) 

where for a mapping / : R d -> M d , V/ = [V/i, . . . , V/ rf ] T . Note that after linearization extra additive 
terms appear as seen in the setting described by equations dU)-©. 



2.2 Message passing 

Assume at time n, the estimate of the localization parameters is 9 n = {^n J }(ij)e£') with 9 % n known to 
node j only. To perform the prediction and update steps in (IT])-® locally at each node a naive approach 
might require each node to access to all localization parameters 8 n and all the different model parameters 
{(C^, d r n , i?^)} n >i. re y . A scheme that requires all this information to be passed at every node would 
be inefficient. It would require a prohibitive amount of communication even for relatively few nodes and 
redundant computations would be performed at the different nodes. The core idea in this paper is to 
avoid this by storing the parameters in 9 n across the network and perform required computations only at 
the nodes where the parameters are stored. The results of these computations are then propagated in the 
network using an efficient message passing scheme. 

Message passing is an iterative procedure with k = 1, . . . ,K iterations for each time n and is steered 
towards the development of a distributed Kalman filter, whose presentation is postponed for the next 
section. In Algorithm Q] we define a recursion of messages which are to be communicated between all pairs 
of neighboring nodes in both directions. Here ne(i) denote the neighbors of node i excluding node i itself. 
At iteration k the computed messages from node i to j are matrix and vector quantities of appropriate 
dimensions and are denoted as rn^, and rfi^ k respectively. The source node is indicated by the first letter 
of the superscript. Note that during the execution of Algorithm Q] time n remains fixed and iteration k 
should not be confused with time n. Clearly we assume that the sensors have the ability to communicate 
much faster than collecting measurements. We proceed with a simple (but key) lemma concerning the 
aggregations of sufficient statistics locally at each node. 



6 



Lemma 1 At time n, let {F^} v( - V be a collection of matrices where F% is known to node v only, and 
consider the task of computing EueV and E^ev F^n a t each node r of a network with a tree topology. 
Using Algorithm 1 and if K is at least as large as the number of edges connecting the two farthest nodes 
in the network, then Y^veV K = K + E m n,K and E uG v F n°n V = E m£ K . 

(The proof, which uses ([5]), is omitted.) An additional advantage here is that if the network is very 
large, in the interest of speed one might be interested in settling with computing the presented sums only 
for a subset of nodes and thus use a smaller K. This also applies when a target traverses the field of view 
of the sensors swiftly and is visible only by few nodes at each time. Finally, a lower value for K is also 
useful when cycles are present in order to avoid summing each more than once, albeit summing only 
over a subset of V. 



3 Distributed Joint Filtering 

For a linear Gaussian system, the joint filter Y\- n ) at node r is a Gaussian distribution with a specific 

mean vector \f n and covariance matrix YT n . The derivation of the Kalman filter to implement pg(x^|Yi :n ) is 
standard upon noting that the measurement model at node r can be written as Y n = C n X^+d n + W n where 
the i-th block of Y n , Y£, satisfies Y* = C % n {X T n + 6 r ' 1 ) + d l n + W^. However, there will be "non-local" steps 

due to the requirement that quantities £ (Q) T ( J R;)" 1 Q, £ (Q) T \R n )~ l Yl and £ (C£) T (i£) _1 Cj;0 r ' i 

iev ieV iev 

be available locally at node r. To solve this problem, we may use Lemma Q] with F l n = (C*) T '(i^)" 1 Q 

and in order to compute ^ (C^) T (-R^) -1 Y^ we will define rn^ k that is an additional message similar to 

iev 

m n,k- 

Recall that b\ , d l n are known local variables that arose due to linearization. Also to aid the development 
of the distributed on-line localization algorithms in Section [H we assume that for the time being the 
localization parameter estimates {# n }n>i are time- varying and known to the relevant nodes they belong. 
For the case where that b l n ,d l n = 0, we summarize the resulting distributed Kalman filter in Algorithm [2 
which is to be implemented at every node of the network. Note that messages (|18p -([20|) are matrix and 
vector valued quantities and require a fixed amount of memory regardless of the number of nodes in the 
network. Also, the same rule for generating and combining messages are implemented at each node. The 
distributed Kalman filter presented here bears a similar structure to the one found in |30| . However, the 
message passing scheme is different and due to the application in mind we have extra terms relevant to 
the localization parameters. 

In the case ^ modifications to Algorithm [2] are as follows: in (|2ip . to the right hand side of 

^n|n-i' ^ ne term b r n should be added and all instances of Y£ should be replaced with Y£ — d r n . Therefore 
the assuming b l n , d l n = does not compromise the generality of the approach. A direct application of this 
modification is the distributed EKF, which is obtained by adding the term (^(/z^^) — V</) n (/i^_ 1 )//^_ 1 
to the right hand side of M^| n _ 1 in (|2ip . and replacing all instances of Y£ with Y£ — V'n^nln-i) ~'~ 
Vipn(fJ! n i n -i)f J ' 1 n i n -i- 111 addition, one needs to replace A n with V0 n (/i^_ 1 ). 

4 Distributed Collaborative Localization 

Following the discussion in Section [2] we will treat the sensor localization problem as a static parameter 
estimation problem for HMMs. The purpose of this section is to develop a fully decentralized implemen- 
tation of popular Maximum Likelihood (ML) techniques for parameter estimation in HMMs. We will 
focus on two on-line ML estimation methods: Recursive Maximum Likelihood (RML) and Expectation- 
Maximization (EM). For the sake of completeness, we have added brief descriptions of these techniques in 
Section [7] of the appendix. 
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Algorithm 2 Distributed Filtering 



1: begin 

2: for n > 1: 

3: Let the localization parameter be 6 n and the set of collected measurements be Y n = {Y^}„ g y • Initialize 
messages { m% nki™nki™nk) anc ^ ( m n fe' fc> fc) ^ or a ^ neighboring nodes S £ as: 



<1 = (Cn) T (K)~ 1 Cn 1 
= (C) T (K)- 1 ^, 



4: for k = 2, . . . , K exchange the messages {mlfk>™l?k'™nk) an< ^ ( m n\>™n\>™n % k) defined below be- 
tween all neighboring nodes € £: 

m ni = (^n) T (K)~ 1 Cn + £ <U ( 18 ) 

p£ne(i)\{j} 

p€ne(i)\{j} 

= + E C 20 ) 

pene(i)\{j} 



5: end for 

6: Update the local filtering densities at each node r £ V: 



Mnln-l — AiMn-lJ ^n|n-l ~~ ^n^n-l^n + Qnj (21) 

k = (s;,^)- 1 + (csfTO- 1 ^ + E m « r ( 22 ) 

igne(r) 

z n = (^n|n-l) Vn|n-1 + (^^(-^n) (23) 

+ e K r -^ r ), 

igne(r) 

ejxK) -1 , ^ = % (24) 



7: end for 
8: end 



S 



The core idea in our distributed ML formulation is to store the parameter 6 n = {9n 3 }(i,j)e£ across 
the network. Each node r will use the available data Y\. n from every node to estimate 9l' 3 , which is the 
component of 9* corresponding to edge (r,j). This can be achieved computing at each node r the ML 
estimate: 

Q r n = arg max logp r 9 (Y l:n ). (25) 

Note that each node maximizes its "local" likelihood function although all the data across the network is 
being used. 

On-line parameter estimation techniques like the RML and on-line EM are suitable for sensor local- 
ization in surveillance applications because we expect a practically indefinite length of observations to 
arrive sequentially. For example, objects will persistently traverse the field of view of these sensors, i.e. 
the departure of old objects would be replenished by the arrival of new ones. A recursive procedure is 
essential to give a quick up-to-date parameter estimate every time a new set of observations is collected 
by the network. This is done by allowing every node r to update the estimate of the parameter along edge 
(r, j), (fn , according to a rule like 

Cii = G£ii(0n,*g, n>l, (26) 

where G^.y is an appropriate function to be defined. Similarly each neighbor j of r will perform a similar 
update along the same edge only this time it will update 9lf . While updating both parameters associated 
to each edge is redundant, it allows a fully decentralized implementation since no other communication is 
needed other than the messages defined in Algorithm [TJ Alternatively one could assign both parameters of 
an edge to just one controlling node. For example in the three node network of Figure [U the parameters 
of edge (1, 2), 9 n ' 2 and 9 n ' , could be assigned to node 2, with the latter having at each time n to update 
9" 2 ,' 1 using an expression like (|26p and then send 9 n ' 2 = —9n l to node 1. 

4.1 Distributed RML 

For distributed RML, each node r updates the parameter of edge (r, j) using 

(27) 

9—9 n 

where 7^ +1 is a step-size that should satisfy Ylmln = 00 an d J2n (in) 2 < °°- 

The gradient in (|27p is w.r.t. O 1 "' 3 . The local joint predicted density pg(x r n \Yi :n -i) at node r was 
defined in ([8]) and is a function of 9 = {# M }(jj)££, and likelihood term is given in ([9]). Also, the gradient 
is evaluated at 9 n = {9 % n}(ij)^e while only On 3 is available locally at node r. The remaining values 9 n are 
stored across the network. All nodes of the network will implement such a local gradient algorithm with 
respect to the parameter associated to its adjacent edge. We note that (|27|) in the present form is not an 
on-line parameter update like (126|) as it requires browsing through the entire history of observations. This 
limitation is removed by defining certain intermediate quantities that facilitate the online evaluation of 
this gradient in the spirit of |18[ [T9] (see in the Appendix for more details). 

The distributed RML implementation for self-localization and tracking is presented in Algorithm [3l 
while the derivation of the algorithm is presented in the Appendix. The intermediate quantities (|28p -([30|) 
take values in R 2 and may be initialized to zero matrices. For the non-linear model, when an EKF 
implementation is used for Algorithm (2j then Algorithm [3] remains the same. 

4.2 Distributed on-line EM 

We begin with a brief description of distributed EM in an off-line context and then present its on-line 
implementation. Given a batch of T observations, let p be the (off-line) iteration index and 9 p = {0p 3 }(i.j)e£ 



7 n+l 



y;j 



+ 7 



n+l 



Vflr.j iog / ^(y„K)p0«|yi :n _i)d< 
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Algorithm 3 Distributed RML 
1: begin 

2: for n > 1: let the current parameter estimate be 9 n . Upon obtaining measurements Y n = {Y^}„ e y 

the following filtering and parameter update steps are to be performed. 
3: Filtering step: Perform steps (3-6) in Algorithm [2j 

4: Parameter update: Each node r G V of the network will update the following quantities for every 
edge (r, j) G 8: 

ii%_ x = A n ^ (28) 

An J = (K)" 1 ^'- (30) 

Upon doing so the localization parameter is updated: 

C'ii = # + 7; +1 [-(^f n _ 1 ) T (^ |n „ 1 )"V;i n _ 1 

+ (i^) T (M;)- 1 4 + m^-m^]. 



5: end for 
6: end 



be the current estimate of 0* after p — 1 distributed EM iterations on the batch of observations Y\-t- Each 
edge controlling node r will execute the following E and M steps to update the estimate of the localization 
parameter for its edge. For iteration p = 1,2,... 

(E step) Q r (9 p ,9) = J logp r e (x r UT , Y 1:T )p r gp (x r 1:T \Y 1:T )dx r 1:T , 
(M step) 9% = argmax Q r (9 p , (9^, 9^)), 

where 9p (rj) = {9^} ee£ \ {rJ y 

To show how the E-step can be computed we write Pg(x^. T , Y\ : t) as, 



T 



Pe( x l:T)Pe( Y l.T\x r i:T) = Yl K-l )P0 ( Y n \ x 



r ) 
nJi 



n=l 



where p r g(Y n \x T n ) was defined in ([9]). Note that Pg p (x\. T \Y\ : T) is a function of 9 p = {9p l }(j j/) g £ (and not 

just Op 3 ) and the #-dependance of Pg(x^. T , Y± : t) arises through the likelihood term only as Pgix^.j.) is 
^-independent. This means that in order to compute the E-step, it is sufficient to maintain the smoothed 
marginals: 

Pe( x n\ Y i:T) oc J p r e (x r 1:T , Y 1:T )dx r l:T y {n} , 

where 1 < n < T and dx^.j,^^ means integration w.r.t. all variables except x r n . For linear Gaussian models 
this smoothed density is also Gaussian, with its mean and covariance denoted by respectively. 

The M-step is solved by setting the derivative of Q r (9 p , (0 r,3 \0 p ^ r '^)) w.r.t. r ' J to zero. The details 
are presented in the Appendix and the main result is: 

V 9 r,j / logpl(Y n \x r n )p r ep (x r n \Y 1:T )dx r n = m?£ K - rh 3 £ K - (m J ^ K ) T fi r n]T , 
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where {j^nKi^nKi^nK)i defined in (fT8|) -(f20j). are propagated with localization parameter 9 p for all 
observations from time 1 to T. Only m^K is a function of Q r > 3 . To perform the M-step, the following 
equation is solved for 6 r ' J 

(£ <KW' j = E«K - «K) T Sn\T ~ <K + <l)- (31) 
n=l n=l 

Note that Q r ' 3 is a function of quantities available locally to node r only. The M-step can also be written 
as the following function: 

Kfjqrj qr,j q r,j v _ / ~r, j \ ~ ( q r,j _ ~r,j\ 
J H°T,1' °T,2> °T,3/ ~~ I °T,2 J l°T,3 ' 

where Sj^, Sjt 3 3 are three summary statistics of the form: 




with s^y, being defined as follows: 

r,j i r v \ ■ j,r .. j,r .. j,r 

s n,3(. x n) r nJ — m n,K ~ m n,K ~r m n,l- 

Note that for this problem s^' J 2 and s^g are state independent. 

An on-line implementation of EM follows by computing recursively running averages for each of the 
three summary statistics, which we will denote as S^\, S^%, $n 3- At each time n these will be used at 

every node r to update r,J using = A(«S^, ^ J 3 ). Note that A is the same function for every 
node. The on-line implementation of distributed EM is found in Algorithm [JJ All the steps are performed 
with quantities available locally at node r using the exchange of messages as detailed in Algorithm [2j 
The derivation of the recursions for S^\, <^n3 are based on (|42[) -([43|) in the Appendix. Here ^ r n is 

a step-size satisfying the same conditions as in RML and 8q can be initialized arbitrarily, e.g. the zero 
vector. Finally, it has been reported in |31] that it is usually beneficial for the first few epochs not to 
perform the M step in (132|) and allow a burn-in period for the running averages of the summary statistics 
to converge. 



5 Numerical Examples 

The performance of the distributed RML and EM algorithms are studied using a Linear Gaussian and a 
non- linear model. For both cases the hidden target is given in ([1]) with V n = BV n , where V n is zero mean 
Gaussian additive noise with variance Q n , and 
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and I is the identity matrix. For the linear model the observations are given by ([2]) with 
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Algorithm 4 Distributed on-line EM 
1: begin 

2: for n > 1: let the current parameter estimate be 9 n . Upon obtaining measurements Y n = {Y^}„gy 

the following filtering and parameter update steps are to be performed. 
3: Filtering step: Perform steps (3-6) in Algorithm [2j Also compute 



4: Parameter update: Each node r £ V of the network will update the following quantities for every 
edge (r, j) £ £: 



-1 




s: 



n,l 



//;;•'//;, • /<;/•'• 



s: 



f,3 
n,2 



ln m n r ,K + C 1 - 7n) S \ 



n-1,2) 



5; 



r t ■ h r ■■ h r 1 ■■ Ji r \ 1 / 1 r\ c' 
7n(™„> - K> + + ( X ~ 7n) 5 , 



n-1,3' 



Upon doing so the localization parameter is updated: 



e 



n+l ~ 



;,l'°n,2'°n,3/- 



(32) 



5: end for 
6: end 
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where a r are constants different for each node and are assigned randomly from the interval [0.75,1.25]. 
For the non-linear model we will use the bearings-only measurement model. In this model at each node 
r, the observation Y£ is: 

Y^ = tan-\Xl{l)/X^)) + W:. 

with W£ A/"(0, 0.35 2 ). For the remaining parameters we set r = 0.01, a x =1 and 9q 3 = for all 
(r, j) E £ . In Figure [2] we show three different sensor networks for which we will perform numerical 
experiments. 

In Figure [3] we present various convergence plots for each of these networks for a y = 0.5. We plot both 
dimensions of the errors 0l — Ol' 3 for three cases: 



• in (a) and (d) we use distributed RML and on-line EM respectively for the network of Figure 2(a) 
and the linear Gaussian model. 

• in (b) and (e) we use distributed RML for the bearings only tracking model and the networks of 
Figures 2(a) and |2(bj] respectively. Local linearization as discussed in Sections 12.11 [31 and 14.11 was 



used to implement the distributed RML algorithm. We remark that we do not apply the online EM 
to problems where the solution to the M-step cannot be expressed analytically as some function A 
of summary statistics. 



• in (c) and (f) we use distributed RML and on-line EM for respectively for the network of Figure 2(c) 
and the linear Gaussian model. In this case we used K = 2. 

All errors converge to zero. Although both methods are theoretically locally optimal when performing 
the simulations we did not observe significant discrepancies in the errors for different initializations. For 
both RML and on-line EM we used for n < 10 3 a constant but small step-size, 7^ = 7 = 4 x 10 -3 and 
0.025 respectively. For the subsequent iterations we set "f n = 7(71 — 10 3 )~ ' 8 . Note that if the step-size 
decreases too quickly in the first time steps, these algorithms might converge too slowly. In the plots 
of Figure [3] one can notice that the distributed RML and EM algorithms require comparable amount of 
time to converge with the RML being usually faster. For example in Figures [3] (a) and (d) we observe 
that RML requires around 1000 iterations to converge whereas on-line EM requires approximately 2000 
iterations. We note that the converge rate also depends on the specific network used, the value of K and 
the simulation parameters. 

To investigate this further we varied K and ^ and recorded the root mean squared error (RMSE) for 
6 n obtained for the network of Figure [2(b) | using 50 independent runs. For the RMSE at time n we will 



USe V 50|£| See£ Em= 



where Qn,m denotes the estimated parameter at epoch n obtained 
from the m-th run. The results are plotted in Figure U] for different cases: 

• in (a) and (b) for ^ = 2 we show the RMSE for K = 2,4, 8, 12. We observe that in every case the 
RMSE keeps reducing as n increases. Both algorithms behave similarly with the RML performing 
better and showing quicker convergence. One expects that observations beyond your near immediate 
neighbors are not necessary to localize adjacent nodes and hence the good performance for small 
values of K. 

• in (b) and (c) we show the RMSE for RML and on-line EM respectively when 2* = 10, 1, 0.5, 0.1. We 
observe that EM seems to be slightly more accurate for lower values of ^ with the reverse holding 
for higher values of the ratio. 

In each run the same step-size was used as before except for RML and 2a = 10, where we had to reduce 
the step size by a factor of 10. 
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(a) 11 node sensor network (b) 44 node sensor network (c) 11 node sensor network with cy- 

cles 

Figure 2: Various sensor networks of different size and topology. 






(a) RML for tree network 



(b) Nonlinear RML for tree network (c) RML for network with cycles 






(d) EM for tree network 



(e) Nonlinear RML for large network (f) EM for network with cycles 



Figure 3: The convergence of the localization parameters' estimates to is demonstrated using appro- 
priate error plots for various sensor networks. Left: Parameter error after each iteration for each edge of 
the medium sensor network of Fig. 2(a) In each subfigure left and right columns show the errors in the 



x- and y- coordinates respectively; (a) is for RML and (d) is for EM. Middle: Same errors when using 



RML for the nonlinear bearings-only observation model; (b) is for medium sized network of Fig. 2(a) and 
(e) for the large network of Fig. |2(b)| Right: Same errors for network with cycles seen in Fig 2(c) ; (c) for 
RML and (f) for EM. 
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(a) RML for different K 



(b) EM for different K 



(c) RML for different ^ 



EM for different 2* 



Figure 4: Comparison of distributed RML and on-line EM. (a) (and (b) resp.): RMSE for RML (and 
on-line EM resp.) against n for K = 2 (□), 4 (o), 8 (o), 12(x). (c) (and (d) resp.): RMSE for RML (and 
on-line EM resp.) for §» = 10 (□), 1 (o), 0.5 (o), 0.1(x). 

6 Conclusion 

We have presented a method to perform collaborative tracking and self-localization. We exploited the 
fact that different nodes collect measurements of a common target. This idea has appeared previously in 
[9j [10], both of which use a Bayesian inference scheme for the localization parameters. We remark that 
our distributed ML methods appear simpler to implement than these Bayesian schemes as the messages 
here are nothing more than the appropriate summary statistics for computing the filtering density and 
performing parameter updates. There is good empirical evidence that the distributed implementations of 
ML proposed in this paper are stable and do seem to settle at reasonably accurate estimates. A theoretical 
investigation of the properties of the schemes would be an interesting but challenging extension. Finally, 
as pointed out by one referee, another interesting extension would be to develop consensus versions of 
Algorithm 1 in the spirit of gossip algorithms in [32j or the aggregation algorithm of |33] which might be 
particularly relevant for networks with cycles, which are dealt with here by using an appropriate value for 
K. 

7 Maximum likelihood parameter estimation 

This section does not pertain to sensor localization specifically but to the general problem of static pa- 
rameter estimation in HMMs using ML. Thus to avoid confusion with the localization problem a different 
font is used the notation. Consider a HMM where {X n } n >i is the hidden state-process and {Y n } n >i is 
the observed process each taking values in taking values in M. dx and W dy respectively. For the transition 
density for {X n } n >i, we have X n+ i|X n = x n ~ f(-|x n ). The observation model, Y n |X n = x„ ~ g$H x n) 
is parametrized by •§ G O (c W 1 " 8 ). The true static parameter generating the sequence of observations is 
i?* and is to be learned from the observed data {Y n } n >i. The ML parameter estimate is the maximizing 
argument of the log-likelihood of the observed data up to time n: i9 n = argmax# 6 @ logp#(Yi :n ). Here 
Pi?(Yi :n ) denotes the joint density of Yi ;n and the subscript makes explicit the value of the parameter used 
to compute this density. 

For a long observation sequence we are interested in a recursive parameter estimation procedure in 
which the data is run through once sequentially. If i? n is the estimate of the model parameter after n 
observations, a recursive method would update the estimate to $ n +i after receiving the new data Y n . For 
example, consider the following update scheme: 

tf n+ i = G n+1 ($ n ,Y n ), n>l. (33) 

where is an appropriate function to be defined. This scheme was originally suggested by |34|. I35| when 

{X n } n >i is not a Markov chain but rather an independent and identically distributed (i.i.d.) sequence. 
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7.1 Recursive Maximum Likelihood (RML) 

To motivate a suitable choice for G n +i($ n , Y n ) for estimating the parameters of a HMM, consider the 
following recursion: 

$n+l = An + ln+1 V logPtf (Y n |Yi :n _i) |^ = ^ n ■ (34) 

where {7«} is the step-size sequence that should satisfy the following constraints: X^n. 7« = 00 ana - 
Xm7n < 00 • O ne possible choice would be j n = n~ a , 0.5 < a < 1. Here Pi?(Y n |Yi :n _i) is the conditional 
density of Y n given Yi :n _i and the subscript makes explicit the value of the parameter used to compute 
this density. Upon receiving Y n , i9 n is updated in the direction of ascent of the conditional density of this 
new observation. The algorithm in the present form is not suitable for online implementation due to the 
need to evaluate the gradient of logp,?(Y n |Yi :n _i) (w.r.t. 1?) at •& = "& n . Doing so would require browsing 
through the entire history of observations. This limitation is removed by defining certain intermediate 
quantities that facilitate the online evaluation of this gradient |18[ PT9] . 

In particular, assume that from the previous iteration of the RML, one has computed p n (x n ) ~ 
P0(xn\Yl:n-l)\#=# n and Pn(*-n) ~ Vp. 1? (x n |Yi :n „i)| 1?= ^ , where (p n ,Pn) are approximations of the pre- 
dicted density and its gradient evaluated at 1? = t? n . The RML is initialized with an arbitrary value for 
$1, Pi(xi) = p 1 ? 1 (xi), which is the prior distribution for Xi and pi(xi) = Vp^(xi)| 1?=1 j , i.e. the gradient 
of this prior which could be zero if it does not depend on Then the online version of (134j) . which is the 
RML procedure of |18l I19| . proceeds as follows. Given the new observation Y n , update the parameter: 

i? n+ i = -d n + 7n+i ( / gtf n (Y n |x n )p n (x n )<ix n ) ( / g,j n (Y n |x n )p n (x n )cZx n + / g^ n (Y n |x n )p n (x n )dx n 



(35) 

where n > 1 and g$'(y|x) = V#g$(y|x)L = ^,. In (f34|) . the desired gradient is the ratio of the terms 
Pi?(Yn|Yi :n _i)|^ =1? and Vp^(Y n [Yi :n _i)|^ =1? . This ratio is approximated in the fraction on the right- 
hand side of (f35|) . After computing (|35|) . one may update (p n ,p n ) to (p n +i,p n +i) for the next RML 
iteration. Specific expressions for this update may be found for example in |141 Section 8.2.1] or [18]. The 
recursive propagation of (p^, p n ) implicitly involves the previous values of the parameter, i.e. and 
hence are only approximations to p ! ?(x n |Yi :n _i)| l?=1?n ^ Vp^(x n |Yi :n _x)|^ = ^ n respectively. It has been 
shown in [18] that the solution of RML converges to the true ML estimator without any loss of efficiency. 
For more details on the convergence of RML for HMMs we refer the reader to |18| . 

7.2 On-line Expectation-Maximization (EM) 

We begin this section with a brief description of Expectation-Maximization (EM) |36| and then present its 
on-line implementation. EM is an iterative off-line algorithm for learning #*, which consists of repeating 
a two step procedure given a batch of T observations. Let p be the (off-line) iteration index. The first 
step, the expectation or E-step, computes 

Q('& p ,'&) = J logp,?(x 1:T , Y 1: T)p0p(xi:T|Yi:T)dXl:T- (36) 

The second step is the maximization or M-step that updates the parameter # p , 

# p+ i = argmax Q(i9 p ,i?) (37) 

Upon the completion of an E and M step, the likelihood surface is ascended, i.e. p# +1 (Yi : r) > p& (Yi-.t) 
|36| . When p$(x± : T, Yi-.t) is in the exponential family, which is the case of linear Gaussian state-space 
models, this procedure can be implemented exactly. Then the E-step is equivalent to computing a summary 
statistic of the form 

Si p = ^ / I \ 's„(x n „i :n ,Y n ) I p^Jx 1:T \Y 1:T )dx 1:T . (38) 



>T P = J, I ^5Z Sn ( X «-l:"> Y ™)^ Ptf P ( X l:T|Yl :T )dx 1:T . 
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where s n : R dx x R d * x R d y ->• M K . In addition, the maximizing argument of Q(i? p ,#) can be characterized 
in this case explicitly through a suitable function A : R K — >• 0, i.e. 

#p = A (sj) . (39) 

Note that in the usual EM setup one has to compute ([380 for every iteration p of the algorithm. 

It is also possible to propose an on-line version of the EM algorithm. This was originally proposed for 
finite state-space and linear Gaussian models in |21[ [371 [20] and for exponential family models in |22[ [3T] . In 
the online implementation of the EM, running averages of the sufficient statistics are computed |20 [ [2~T [ [22] . 
Let {$m}i<m<n be the sequence of parameter estimates of the online EM algorithm computed sequentially 
based on Yi :n _i. When Y n is received, we compute 

= 7n f S n (Xn-hn) P# 1: n y^n— l:n 

|Yl :n )dx n _l : „ 

+ (1 - In) Em=\( EI (1 - li))lm /Sra (Xm-l:m)W 1: „ (x m _i :m |Yi :n )dx m _i :m , 
i=m+l 

where the subscript $i :n on p 1 j 1 . n (xi : T|Yi :n ) indicates that the posterior density is being computed sequen- 
tially using the parameter # m at time m < n. The step sizes {7n} n >i need to satisfy YJ n 7„ = oo and 
7n < 00 as m the RML case. For the M-step one uses the same maximization step ([39]) used in the 
batch version 

tf n+ i = A(<S n ). (41) 
The recursive calculation of S n can be achieved by setting V\ (xq) =0 and computing 



V n (x n ) = J {j n S n (x n „!,X n ) + (1 - 7n) K-l (Xn-l)} 

:n— 1) X« ) dx n _! (42) 
K (x n )^ 1:n (x n |Y 1:n )(ix n . (43) 



and 



For finite state-space and linear Gaussian models, all the quantities appearing in this algorithm can be 
calculated exactly |20l |2T[ l3T] . 



8 Distributed RML derivation 

Let 9 n = {On' 1 } (ij)ee be the estimate of the true parameter 8* given the available data ii :n _i. Consider 
an arbitrary node r and assume it controls edge (r, j). At time n, we assume the following quantities are 
available: = Vgr.j p<n-l \ $—$ > ^n-l\e—e i^n-l)- The ^ rs ^ °f these quantities is the derivative of the 

conditional mean of the hidden state at node r given Yi :n _i, i.e. V e r, 3 f x 7 n _ 1 p r e (x 1 n _ l \Yi- n -i)dx 7 n _ 1 \ e=d . 
This quantity is a function of the localization parameter 9 n . is the variance of the distribution 

PeC^n-ll^n-l) L_g an d is independent of the localization parameter. The log-likelihood in ([27]) evaluates 
to: 

\0gp r e {Yn\Yl:n-l) = "^E^" " K^K ~ ^) 

iev 

~ ^n|n-l T (S„|n-l)~V„|n-l + 5«) T (K) _1 < + 

where all independent terms have been lumped together in the term 'const'. (Refer to Algorithm [2] for 
the definition of the quantities in this expression.) Differentiating this expression w.r.t. # r ' J yields 

v^iogp^iYi^) = -(v^/x; |n _ 1 ) T (E; |ri _ 1 )-V^_ 1 

+ (v^o t (m;)- i 4 + Y^(y^P'Y(cti T &rHK - cie r n. 

iev 
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(I27|) requires Vgrj log pg(Y n \Yi :n _i) to be evaluated at 9 = 9 n . Using the equations (|2ip -([24")) and the 

assumed knowledge of (fi^-ii Hn-ilg— q ; ^n-l) we can evaluate the derivatives on the right-hand side of 
this expression: 



,\'<J — V7 . ,, r 



An^-V (44) 



= V^K|, =0n = (M;)- 1 ^. (46) 

Using property ([5]) we note that for the set of vertices i for which the path from r to i includes edge 
(r, j), V ' Qr,j9 r '' 1 = I (the identity matrix) whereas for the rest Vgr,j#™ = 0. For all the nodes i for which 
V0r,j# r ' 1 = I, let them form a sub tree (V' r j,£' r j) branching out from node j away from node r. Then the 
last sum in the expression for Vgrj log ^(l^Yi^-i)^^ evaluates to, 



m n,K ~ m n,K> 



where messages (m^K'^ii k) were defined in Algorithms [2 Similarly, we can write the sum in the 
expression for Zn as m^ r K (again refer to Algorithms [2]) to obtain 



^ = (EWXfn- 1 -<*:- (47) 



To conclude, the approximations to (/in"' = Vgrj/x^| 6 , =6lji+i , ^ r n \g = g n+1 ,^n) f° r the subsequent RML itcr- 



{ r,3 

ation, i.e. (|27f> at time n = n + 1, are given by 



while (^n\e=e n+1 >^n) are gi ven by (|21|) - f|24=[) . The approximation to ^ e r ^ l JL n\e=9 n+1 f°h° ws from differ- 
entiating (f24"|) . (Jin 3 , t J 'n\e=e n+1 ) are only approximations because they are computed using the previous 
values of the parameters, i.e. 9\ :n . 



9 Distributed EM derivation 

For the off-line EM approach, once a batch of T observations have been obtained, each node r of the 
network that controls an edge will execute the following E and M step iteration n, 

Q r (9 p ,9) = J logp^^ :T ,F 1:T )p^(^ :T |y 1:T )^ :T; 
6% = arg max Q r (9 p , (9^ , {9 e ,e € £\(r,j)})), 

where it is assumed that node r controls edge (r,j). The quantity p r e (x\-t\Yi\t) is the joint distribution 
of the hidden states at node r given all the observations of the network from time 1 to T and is given up 
to a proportionality constant, 

T 

Pe P ( x v.T)Pe p O r i-.T\x r 1 . T ) = \\ /nC^K-i^P^K), 

n=l 

where pg p (Y n \x n ) was defined in ([9]). Note that Pq p (x\. t , Y\.t) (and hence P t q (x\. t \Y\ : t)) is a function of 
9 P = {9p l }(i t i')££ and not just 9 P :> '. Also, the ^-dependence of Pq(xj, t , Yi : t) arises through the likelihood 
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term only as is 9 independent. Note that 

E log^ocK + o r ' v ) = E < ~ \ E( y n - c^) t (K)-Hy: - cy^) 



2 

ugv vev vev 



vev 



Era T (Kr^ 



vev 



where c" n is a constant independent of 0. Taking the expectation w.r.t. p r Q n {x r n \Yi-.T) gives 

J logpl(Y n \xM p (x r n \Y 1:T )dx r n = -\ E - ^^) T (K)- 1 « - C^'")] 

" (^|T) T E( C n) T «)" 1C n^ + Const 



vev 



where all terms independent of 9 r ' J have been lumped together as 'const' and [i r n \ T is the mean of x r n under 
p T e {x r n \Y\-T)- Taking the gradient w.r.t. r > 3 and following the steps in the derivation of the distributed 
RML we obtain 



V er j J log^(F n |a;;)^ p «|yi : T)d< = m>£ K - m>£ K - (m£ K ) ff n \ T 

where ijn d r ^ K ^rn^ K ^w?^ K ) is defined in (fT8|) -(f20]). Only rh 3 ' r K is a function of 9 r ' J . Now to perform the 
M-step, we solve 



T 



E < K o rJ = E < K - «k) t ^t - E 

\n=l J n=l \ j'Gne(jr)\{r} 



Note that # rj can be recovered by standard linear algebra and so far 9 r ' 3 is solved by quantities available 
locally to node r and j. One can use the fact that ™n K-i = ^nK ~ '"nl ^° so that the 

j'Gne(j)\{r} 

M-step can be performed with quantities available locally to node r only. Recall that 'Y2m=i f^n\T — 

f (Yln=i x n) Po p ( x i-T\Yi--T)dx r 1 . T . This implies directly that three summary statistics are needed for node 

r to update # rj . These should be defined using: 

r,j / r y \ / j,r \T r r,j / r \r \ j,r r,j '• / r \r \ ■ j,r ■■ j,r . ■■ j,r 

S n,l\ X m in) — \ m n ^K) X ni S n,2V x ni r n) — m nj K' S n,3\ X m 1 n) — m n ,K ~ m n,K "+~ m n,l> 

where s r nl , s r n3 are each functions of x r n and Y n via fJ^i T and m^'^ — rh'lfx + respectively. The 

summary statistics can be written in the form of (|38p as follows: 

T 



S T^ = \\ EKV) T ^ ) PlM-.T\Y,:T)dxl. T 



\n=l 

T T 

S T,2 = fl^ m n,K' S T% = f 2^ [ m n,K ~ m n,K + m n,lJ > 
n=l n=l 

and the M-step function becomes A(si, S2, S3) = s^ 1 (S3 — si) , where si, S2, S3 correspond to each of the 
three summary statistics. Note that A is the same function for every node. 

We will now proceed to the on-line implementation. Let at time n the estimate of the localization 
parameter be 9 n . Following the description of Section [7^21 for every r € V and (r,j) € £, let S^\, ^ J 2 , 
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S^z be the running averages (w.r.t n) for Sj,^ , an d 5y J 3 " respectively. The recursions for 

are trivial: 

S ni = ln m n r ,K + i 1 ~ 7n) 5 n,3 = 7n(^' K ~ ^n,K + ™n' l) + " 7n) 5 n-l,3> 

where {7^} n>1 needs to satisfy ^n>i 7ri = 00 an d Sn>i (7«) 2 < 00 ■ For 5^, we will use (I4"2l) - (l4"3j) . We 
first set V^ 7 (xq) = and define the recursion 

C J 'K) = 7;(<V) T < + (l-7;) / K^i(^i)p5u»K-i| ^1.^)^1- ( 4 8) 

Using standard manipulations with Gaussians we can derive that Pg Vn (2^-1 1 ^4:n-i, ^n) i s itself a Gaussian 
density with mean and variance denoted by /2^(x n ),£^ respectively, where 

= (^n-1 + ^I^n^n) 1 ^ni x n) = ^ f (E^_ x ) + /^nl • 

It is then evident that ()48|) becomes (a;^) = H^ 3 x r n + with: 

^ = 7;(<x) t + (i-7;)^ 1 (s;)" 1 ^q-i, 
# = (i-7D(<'ii(s;)" 1 K_ 1 ) _1 ^-i+^ii) ; 

where Hq J = and /iq J = 0. Finally, the recursive calculation of S n '^ is achieved by computing 

Again all the steps are performed locally at node r, which can update parameter O 7 "' 1 using 0^+1 = 
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